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Complexity measures are essential to understand complex systems and there are numerous defi- 
nitions to analyze one-dimensional data. However, extensions of these approaches to two or higher- 
dimensional data, such as images, are much less common. Here, we reduce this gap by applying the 
ideas of the permutation entropy combined with a relative entropic index. We build up a numerical 
procedure that can be easily implemented to evaluate the complexity of two or higher-dimensional 
patterns. We work out this method in different scenarios where numerical experiments and empiri- 
cal data were taken into account. Specifically, we have applied the method to i) fractal landscapes 
generated numerically where we compare our measures with the Hurst exponent; ii) liquid crystal 
textures where nematic- isotropic- nematic phase transitions were properly identified; Hi) 12 charac- 
teristic textures of liquid crystals where the different values show that the method can distinguish 
different phases; iv) and Ising surfaces where our method identified the critical temperature and 
also proved to be stable. 



I. INTRODUCTION 

Investigations related to the so called complex systems are widely spread among different scientific communities, 
ranging from physics and biology to economy and psychology. A considerable part of these works deals with empirical 
data aiming to extract patterns, regularities or laws that rule the dynamics of the system. In this direction, the 
concept of complexity measures often emerges. Complexity measures can compare empirical data such as time 
series and classify them in somewhere between regular, chaotic or random [1], while other complexity measures can 
differentiate between degrees of correlations [2]. Examples of these measures include algorithmic complexity [H, 
entropies [4], relative entropies [5[, fractal dimensions [6], and Lyapunov exponents These seminal works are 
still motivating new definitions, and today there are numerous definitions of complexity, which have been successful 
applied to different areas such as medicine 0, Q, ecology [10-13], astrophysics [14-16], and music [17, 18]. 

It is surprising that this large number of complexity measures is mainly focused on one-dimensional data, while 
much less attention has been paid to two and higher-dimensional structures such as images. Naturally, there are few 
exceptions such as the work of Grassberger [l9[ and more recent Refs. (2Ql-[22j. though some of the authors of these 
papers agree that a higher-dimensional approach still represents an open and subtle problem. Furthermore, as it 
was stated by Bandt and Pompe [23], most of the complexity measures depend on specific algorithms or recipes for 
processing the data which may also depend on tuning parameters. As a direct consequence, there are huge difficulties 
for reproducing previous results without the knowledge of details of the methods. 

Bandt and Pompe not only raised this problem, but they also proposed an alternative method that tries to overcome 
the previous problems, introducing what they call permutation entropy — a natural con iplexity measure for time series. 
There are many recent applications of this new technique that confirm its usefulness [2^-31]. In particular. Rosso et 
al. [1] have successful applied the Bandt and Pompe ideas together with a relative entropic measure [32] to differentiate 
chaotic time series from stochastic ones. They have constructed a diagram, which was first proposed by Lopez-Ruiz et 
al. [33], (called as complexity-entropy causality plane) by plotting the relative entropic measure versus the permutation 
entropy. Intriguingly, chaotic and stochastic series are located in different regions of this representation space. 

Here, we show that the complexity-entropy causality plane can be extended for higher-dimensional patterns. We 
apply this new approach in different scenarios related to two-dimensional structures and the results indicate that the 
method is very promising for distinguishing between two-dimensional patterns. The following sections are organized as 
follows. Section II is devoted to review briefly the properties of the permutation information-theory-derived quantifiers 
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and the complexity-entropy causality plane, and also to define an appropriate way to generalize these definitions to 
higher-dimensional data. In Section III, we work out several applications based on numerical and empirical data. 
Section IV presents a summary of our results. 



II. METHODS 



The ingenious idea of Bandt and Pompe [23] was to define a measure that may be easily applied to any type of time 
series. The method lies on associating symbolic sequences to the segments of the time series based on the existence 
of local order, and next, by using probability distribution associated to these symbols, to estimate the complexity 
quantifier. For purpose of definition, let us consider a time series {xt}t=i,...,n composed by n elements and also 
d-dimensional vectors {d > 1) defined by 

where 5 = d + 1, . . . , n. Next, for all the {n — d-\-l) vectors, we evaluate the permutations tt = (ro, ri, . . . , rd-i) of 
(0, 1, . . . , d — 1) defined by Xs-^d-i ^ Xs-rd-2 ^ ' ' ' ^ ^s-ri ^ ^s-ro- The d \ possible permutations of tt will be the 
accessible states of the system, and for each state we estimate the ordinal pattern probability given by 

#{s\s < n — d -\- 1; {s) has type tt} 



p{7r) 



1 



where the symbol # stands for the number of occurrences of the permutation tt. Now, we can apply the ordinal 
patterns probability distribution, P = {p(7r)}, to estimate a complexity measure based on some entropic formulation. 

Before advancing, we note that the previous method may be extended to higher-dimensional data structures such as 
images. In order to do this, we consider that the system is now represented by a two-dimensional array 
of size Ux X riy. In analogy to the vector (5), we define dx x dy matrices {dx^dy > 1) given by 
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where = d^^, + 1, . . . , n^^ and Sy = dy^ dy -\- 1^ . . . ^Uy. Next, for all these {ux — dx -\- — dy ^) matrices, 

we evaluate the permutations tt = [(ro,iXo), (n,^o), • • • , (rrf^-i,^io), • • • , (ro,iXd^-i), {ri,Udy-i), . . . , (r^^-i, ^d^-i)] of 

ro 1 d d -1) defined bv 7/^^"''^^-' < ./^"""^^-^ < • • • < o/^'^'^'y-^ < ./y~''^y-^ < < .,^^-^0 < ./^-^o < 
l^u, i, . . . , ax ay ij aennea oy yg^_^^^_^ 1^ ysx-rd-2 - - ^Sx-ri ^ ysx-ro ^ ^ ysx-rdx-i - ysx-rd-2 - 

• • • ^ Vsl-ri — Vsl-ro ' The system can now access {dx dy)\ states for which we calculate the probability distribution 
P = {p{7t)} through the relative frequencies given by 



P{7t) 



#{{sx,Sy)\sx <nx - dx and Sy < Uy - dy ^ 1; jsx.Sy) has type tt} 
{rix - dx ^ !)(% - dy ^1) 



For easier understanding, we illustrate this procedure for a small array in Fig. [H 

Naturally, the order procedure that defines the permutation tt is no longer unique as in the one-dimensional case. 
For instance, instead of ordering the elements of {sx, Sy) row- by-row, we could also order column- by-column. However, 
these other definitions will only change the "name" of the states in such a way that the set P = {p{7t)} will remain 
unchanged. Thus, there is no lost of generalization in assuming a given order recipe for defining tt. 

We note that this procedure is straightforward generalized to accomplish higher-dimensional structures (e.g., the 
volumetric brain images obtained via functional magnetic resonance imaging), and that it recovers the one-dimensional 
case by setting Uy = 1 and dy = 1. Here, for simplicity, we focus our analysis on two-dimensional structures. 

The parameters dx and dy (known as embedding dimensions) play an important role in the estimation of the 
permutation probability distribution P, since they determine the number of accessible states. In the one-dimensional 
case, it is usual to choose d! <C n in order to obtain reliable statistics in the one-dimensional case (for practical 
purposes, Bandt and Pompe recommend d = 3, . . . , 7 [23]). For the two-dimensional case a similar relationship must 
hold, i.e., {dxdy)\ <C rixTiy. 
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FIG. 1: Schematic representation of the construction of the accessible states. In this example we have a 3 x 3 array 
(left panel) and we choose the embedding dimensions dx — 2 and dy — 2. In the right panel we illustrate the construction 
of the states. We first obtain the sub-matrix corresponding to = 2 and = 2 that have as elements (1,2,8,3) and, 
after sorting, this sub-matrix leads to the state "0132". We thus move to next sub-matrix = 3 and Sy — 2 which have 
the elements (2, 1,3,4) and that, after sorting, leads to the state "1023". The last two remaining matrices lead to the states 
"1230" and "0132". Finally, we estimate the probabilities p(7r), that are, p("0132") = 2/4 = 0.5, p("1023") = 1/4 = 0.25 and 
p("1230") = 1/4 = 0.25 which are then used in the equations ([T|) and ([2|), leading to H ^ 0.33 and C ^ 0.27. 



To go further, we need to rewrite the entropic measures used in Refs. [7, '235. The first one is called normalized 
permutation entropy ^] and it is obtained by applying the Shannon's entropy to the probabilities P = {p(7r)}, i.e., 

H[P] = , (1) 

where S[P] = — ^p{7t) logp(7r) and 5'max = ^og[{dxdy)\]. The value of ^max is obtained by considering all the {dxdy)\ 
accessible states to be equiprobable, i.e., P = = l/{dxdy)\. By definition, < H[P] < 1, where the upper bound 
occurs for a completely random array. We expect H[P] < 1 for arrays that exhibit some kind of correlated dynamics. 
The other measure pL,] is defined by 

C[P]=Q[P,Pe]H[P], (2) 

where Q[P, Pe] is a relative entropic metric between the empirical ordinal probability P = {p{7r)} and the equiprobable 
state Pe = l/(dxdy)\. The quantity Q[P, Pe] is known as disequilibrium and it is defined in terms of the Jensen- Shannon 
divergence [3^ (or also in terms of a symmetrized Kullback-Leibler divergence [35]) and can be written as 

Q[P, P.] = S[iP + Pe)/2]- S[P]/2-S[P.]/2 

where 

Qm.. = -\ 1^^^]^ log[(d,d,)! + 1] - 21og[2(d,d,)!] + log[(d,d,)!]| 

is the maximum possible value of Q[P^Pe], obtained when one of the components of P is equal to one and all the 
other vanish. 

The disequilibrium C quantifies the degree of correlational structures providing important additional information 
that may not be carried only by the permutation entropy. In addition, for a given H[P] value there exists a range 
of possible values for C[P] [36]. This is the main reason why Rosso et al. [i] proposed to employ a diagram of C[P] 
versus H[P] as a diagnostic tool, building up the complexity-entropy causality plane. 
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III. APPLICATIONS 

In the following, we will calculate the diagram of C[P] versus H[P] to measure the complexity and to distinguish 
among different two-dimensional patterns. 

A. Fractal Surfaces 

We generate fractal surfaces through the random midpoint displacement algorithm [s^ • This algorithm starts with 
a square. For each vertex, we assign a random value representing the surface height. Next, we add a new point 
located at the center of the initial square. We set the height of this point equal to the average height of the previous 
four vertex plus a Gaussian random number with zero mean and standard-deviation 5i. We also add four points 
located at the middle segments which connects each initial vertex. For these four points, the heights are equal to the 
average value between the two closest vertex and the middle point plus a Gaussian random number with zero mean 
and standard-deviation 5i. Now, we imagine that these 9 points represent four new squares and, for each one, we 
apply the previous procedure using (^2. By repeating this process k times and using 5k = 2~~, we should obtain 
a square surface of side 2^ + 1 with fractal properties. Here, h is the Hurst exponent and D = ?> — h \s the surface 
fractal dimension. Figure [2] shows several surfaces generated through this procedure for different values of h. 



h ^ 02 



h = 0.3 




I 



i 







10 20 30 40 50 60 10 20 30 40 50 60 10 20 30 40 50 60 



FIG. 2: Examples of fractal surfaces obtained through the random midpoint displacement method. These are 
65 X 65 surfaces (/c = 6) for different values of the Hurst exponent h. For easier visualization, we have scaled the height of the 
surfaces in order to stay between and 1. We note that for small values of h the surfaces display an alternation of peaks and 
valleys (anti-persistent behavior) much more frequent than those one obtained for larger values of h. For larger values of /i, the 
surfaces are smoother reflecting the persistent behavior induced by the value of h > 0.5. 
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We apply our method for these surfaces aiming to verify how the permutation quantifiers H and C change with 
the Hurst exponent /i, as it is shown in Fig. [3l In these 3d plots, we show the localization in the causality plane 
obtained for different values of h evaluated from 1025 x 1025 surfaces {k = 10). In Fig.[3K, we use dx = 2 and dy = 3 
(circles), and dx = 3 and dy = 2 (squares) as embedding dimensions. Note that the values of H and C are practically 
invariant under the rotation dx dy and dy dx. This invariance is related to the fact that in these fractal surfaces 
there is not preferential direction. In Fig. [Sd, we employ dx = 3 and dy = 3. We note basically the same dependence 
but a different range for H and C, since this change increases the number of accessible states. These results show 
that our method properly differentiates fractal surfaces concerning the Hurst exponent. Moreover, we investigate the 
robustness of the permutation quantifiers under several realizations of the random midpoint displacement algorithm 
and the results show that both indexes are very stable. For example, the standard-deviation in the values of H and 
C are usually smaller than 10~^ when considering k = 10. 




FIG. 3: Dependence of the complexity-entropy causality plane on Hurst exponent h. We have employed fractal 
surfaces of size 1025 x 1025 {k = 10). In (a) we plot C and H versus h for the embedding dimensions dx = 2 and dy = 3 
(circles) and also for dx = 3 and dy — 2 (squares). We note the invariance of the index against the rotation dx dy and 
dy ^ dx. In (b) we plot the diagram for dx = dy = 3. We observe changes in the scale of C and H caused by the increasing 
number of states. In both cases, as h increases the complexity C also increases while the permutation entropy H decreases. 
This behavior reflects the differences in the roughness shown in Fig.O For values oih < 0.5 the surface is anti-persistent which 
generates a flatter distribution for the values of P = {p{^)} leading to values of C and H closer to the aleatory limit (C ^ 
and H ^ 1). For values of h > 0.5 there is a persistent behavior in the surfaces heights which generates a more intricate 
distribution of P = {pM} and, consequently, values of H and C that are closer to the middle of the causality plane (region of 
higher complexity). 



B. Liquid Crystal Textures 

Another interesting application is related to different patterns that a thin film of a liquid crystal exhibits. These 
textures are obtained by observing a thin sample of liquid crystal placed between two crossed polarizers in a microscope. 
The textures give useful information about the macroscopic structure of the liquid crystal. For instance, different 
phases have different typical textures, and by tracking their evolution one can properly identify the phase transition. 

We first study a lyotropic liquid crystal under isotropic-nematic-isotropic phase transition. Figure 2] shows three 
snapshots of the texture at different temperatures. In this case, we clearly note the differences in the textures. The 
leftmost and rightmost textures are at the isotropic phase while the middle one is at the nematic phase. We observe 
that the pattern is very complex for the nematic phase, while for the isotropic one it is basically random. 

We calculate H and C as a function of the temperature for different values of the embedding dimensions, as it is 
shown in Fig. [5l In these plots, the different shaded regions represent the different liquid crystal phases. We note 
that the phase transitions are successful identified independently of dx and dy. However, Fig[5t and[5li show a slight 
different dependence of H and C versus the temperature when considering dx = 2 and dy = 3 or dx = 3 and dy = 2. 
Because the liquid crystal sample is placed in elongated capillary tube, there is a surface effect that act on the liquid 
crystal molecules. This effect is usually amplified at the phase transition and it is also the reason for differences 
between the embedding dimensions. 
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FIG. 4: Characteristic textures of a lyotropic liquid crystal at different temperatures and phases. The lyotropic 
system used here is a mixture of potassium laurate 27.00%), decanol 6.24%) and deuterium oxide 66.76%) — suitable 
concentrations in order to get a isotropic nematic isotropic phase sequence These images were constructed by 

observing the optical microscopy of a flat capillary which contains the mixture at different temperatures. Here, we have used 
the average value of the pixels of the three layers (RGB) of the original image and a rescaled temperature. 
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FIG. 5: Dependence of the entropic indexes on the temperature of a lyotropic liquid crystal. We plot H versus 
the temperature in (a) and C versus the temperature in (b), where we employ dx — dy — 2. Figures (c) and (d) present the 
results for = 2 and = 3, and also for = 3 and dy — 2. The different shaded areas represent the different liquid crystal 
phases. Note that the phase transitions are properly identified in all cases. Due to the asymmetry of the elongated capillary 
tube where the liquid crystal sample is placed, R and C present slight differences under the rotation dx dy and dy ^ dx. 
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In this particular phase transition, the difference between the textures are large enough that it can be identified 
just by visual inspection. However, this is not the usual case and many phase transitions are very difficult to identify. 
In this context, an interesting question is whether our method can help to distinguish different phases. To address 
this question, we evaluate H and C for twelve characteristic textures of different liquid crystals. We download these 
textures from the webpage of the Liquid Crystal Institute at Kent State University [39] and Fig. [6] shows the value of 
H and C for each texture in the causality plane. The results allow to conclude that the method ranks the textures 
in a kind of complexity order where each characteristic texture occupies a different place in this representation space. 
Moreover, the different values of H and C indicate that the permutation quantifiers can also identify smooth phase 
transitions. 

Naturally, the location of each texture in the causality plane should be related to physical properties of the liquid 
crystals. A better understanding of the relation between the permutation quantifiers and these physical attributes 
may deserves a more careful investigation since some properties of Hquid crystals such as the order parameter can 
be quite hard to empirically measure. In this context, the existence of a clear relation between, for example, the 
order parameter and H or C will be experimentally handy. Here, we just have the pictures of the textures in such a 
way that is very hard to point out these relationships. However, a visual inspection of Fig. [6] suggests that some of 
the more ordered phases, such as the blue phase (this phase display a cubic structure of defects), are located in the 
central part of the causality plane (region of higher complexity) , while other textures which present a large number 
of non-ordered defects, such as the Smectic B and C, are positioned closer to the aleatory limit (C ^ and H ^ 1). 
Thus, it seems that the permutation quantifiers are capturing in somehow the competition between the orientational 
order of the phase and, also, the number of defects present in the textures. 
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FIG. 6: Complexity-entropy causality plane evaluated for several liquid crystal textures |39l |. Here, we have used 
the averaged pixel values of the three layers (RGB) of the original image and dx = 2 and dy = 3. The image sizes are about 
270 X 200 pixels. We note that each texture has a unique position in the causality plane which indicates that the permutation 
quantifiers are capable of differentiate not only transitions involving the isotropic phase, but also smoother phase transitions. 
We further observe that some high ordered phase such as the blue phase are located at the central part of the causality plane 
(region of higher complexity), while other phases which present a large number of defects such as the Smectic B and C are 
closer to the aleatory limit (C ^ and H ^ 1). 
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C. Ising Surfaces 



As a last application, we study the permutation measures H and C applied to Ising surfaces [40|,|41|. These surfaces 
are obtained by accumulating the lattice spin values G { — 1,1} of the Ising model defined by the Hamiltonian 



E 



(JiCFn 



(4) 



where the sum is over all the pairs of first neighbor sites in the lattice. We numerically solve this spin-1/2 Ising model 
on a 1/ X L lattice using the Monte Carlo method with periodic boundary conditions. By using the spin values, we 
define the surface height for each lattice site i as 



(5) 



where t represents the number of Monte Carlo steps. In Fig. [71 we show three surfaces obtained though this procedure 
for different values of the reduced temperature T/Tc, where Tc = 2/ln(l + a/2) is the critical temperature of the 
model. We note the complex pattern exhibited by the surface for T/Tc = 1, and the almost random patterns for 
T/Tc > 1 and T/Tc < 1- 

We first investigate the dependence of H and C on the reduced temperature T/Tc after a large number of Monte 
Carlo steps (10^) and for L = 500. Figures [8^1 and[8]3 show H and C for dx = 2 and dy = 3, and for the rotation 
dx dy and dy dx. We note that, at the critical temperature, both indexes display a sharp peak and that they are 
invariant under the rotation. Moreover, Fig. [8t presents a 3d visualization of the phase transition for dx = dy = 3. 
This higher-dimensional representation can be useful when investigating more complex phase transitions, since a 
greater number of degrees of freedom allows the critical point to be more visible. 

We further study the temporal evolution of H and C for different reduced temperatures, as it is shown in Fig. [9l 
The initial values of the spins were chosen equal to 1 and, as we can see, the values for H and C are different just 
after one Monte Carlo step. For T ^ Tc^ the value of H increases over time and around t ~ 10^ it reaches a plateau. 
For T = Tc^ the value of H increases up to a maximum value around t ^ 10^ and then starts to approach a lower 
plateau value. A striking behavior is observed for C, where for all temperatures the complexity displays a maximum 
value before it begins to approach a plateau value. It is worth noting that both quantifiers are very stable after '^10^ 
Monte Carlo steps. 
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FIG. 7: Examples of Ising surfaces for three different temperatures. These surfaces were obtained after 10^ Monte 
Carlo steps for three different temperatures: below Tc, at Tc and above Tc. In these plots, the height values were scaled to stay 
between and 1. We note that for temperatures higher or lower than Tc, the surfaces exhibit an almost random pattern. For 
values of the temperature closer to Tc the surfaces exhibit a more complex pattern, reflecting the long-range correlations that 
appear among the spin sites during the phase transition. 
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FIG. 9: Dependence of the entropic indexes on the number of Monte Carlo steps. Here, t denotes the number of 
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IV. SUMMARY 

We have proposed a generalization of the complexity-entropy causality plane to higher- dimensional patterns. We 
applied this approach to fractal surfaces, liquid crystal textures and Ising surfaces. It was shown that the indexes 
H and C performed very well for distinguishing between the different roughness of the fractal surfaces. The indexes 
properly identified the phase transitions of a lyotropic liquid crystal and sorted diff'erent characteristic textures in 
a kind of complexity order. Finally, concerning the Ising surfaces, the indexes not only had identified the critical 
temperature, but also proved to be stable after ~ 10^ Monte Carlo steps. The method also has a very fast and simple 
numerical evaluation. Taking into account all these findings, we are very optimist that our method can reduce the 
gap between one- dimensional complexity measures and the higher- dimensional ones. 
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